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ABSTRACT 

A multidomain Chebyshev spectral collocation method for solving hyperbolic 
partial differential equations has been developed. Though spectral methods 
are global methods, an attractive idea is to break a computational domain into 
several subdomains, and a way to handle the interfaces is described. The 
multidomain approach offers advantages over the use of a single Chebyshev 
grid. It allows complex geometries to be covered, and local refinement can be 
used to resolve important features. For steady-state problems it reduces the 
stiffness associated with the use of explicit time integration as a relaxation 
scheme. Furthermore, the proposed method remains spectrally accurate. 
Results showing performance of the method on one- dimensional linear models 
and one- and two-dimensional nonlinear gas-dynamics problems are presented. 
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1. INTRODUCTION 


In this paper we address the problem of efficiently computing Chebyshev 
spectral collocation approximations to quasilinear hyperbolic systems of the 
form 

Q + A(Q)Q + B(Q)Q - 0 x,y DCR 2 , t > 0 (1) 

l x y 

with appropriate boundary and initial conditions. Here, Q is an m-vector 
and A and B are mxm matrices. This system is hyperbolic if for any 
constants kj and k£ the matrix T ■ k^ A + k£ B has only real eigenvalues 
and there exists a similarity transformation matrix, P, such that PTP * ** A 
is a real diagonal matrix. 

In particular, we are interested in the solution of the Euler equations of 
gas dynamics which form a system of this type. The use of the nonconservation 
form is justified for problems in which shocks are fitted and in this 
situation spectral methods work well [1]. Problems of the type presented in 
Ref. [1] provide the motivation for what follows. 

The typical Chebyshev spectral collocation procedure for the solution of 
the system (1) is described in several reviews such as those of Gottlieb, 
Hussaini, and Orszag 12], and Hussaini, Salas, and Zang (3], First, the 
domain of interest is mapped onto the square D' «= [-1 ,l]x[-l ,1] and an 
(N+M) x (M+l) point mesh is generated with the collocation points defined by 


x^ «= - cos(iir/N) i = 0,1,»*»,N 


- cosCjit/M) j *= 0 , 1 , • • • ,M 


( 2 ) 
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Mesh point values of Q, designated by Qfj » are associated with each of the 
collocation points (x^,yj). A global Chebyshev interpolant of order N in 
the x direction and order M in the y direction is then put through the 
mesh point values 

N,M 

Q (*,y> - l a nm T (x)T <y). (3) 
n,m=0 

Approximations to the derivatives at the collocation points are computed 
by differentiating the interpolant and evaluating the resulting polynomial at 
the collocation points. The computation of the derivatives can be 
accomplished in one of two ways (see Gottlieb, et al., [2]): The first is to 
take advantage of the fact that the sums for both the interpolant and its 
derivative reduce to cosine sums at the chosen collocation points. For 
example 


dQ N,M N,M 

dT - l . - l . b.jr.C’OTjy) 


dx “ _ nm ' n 

n,ra=0 


where 


n,m=0 


nm n 


m 


(4) 


Nm 


- o. 


b M . = 2Na 

N-l ,m nm 


(5) 


and 


C n b nm = b n + 2,m + 2(n + 1)a n + l,m for 0 < n < N " 2 - 


The constant c n is defined as c n = 2 for n = 0,N and c n = 1 
otherwise. The advantage of this form is that a fast cosine transform 
compute the derivatives along each y line in 0(N log N) operations. 


can 



-3- 


The other approach to computing the derivatives is to write the 
differentiation operation as the product of a differentiation matrix and the a 
vector of the Q^j's. For example, along each y line the x derivative is 

(^E)J - KQ^ «> 

T 

where (Qp)j “ [Qq j Qj j ••• j] and the elements of the matrix D are 
defined in Gottlieb et al., [2], The amount of work with this procedure is 
of O(N^). What one loses in efficiency one gains as flexibility in the 
number of mesh points that can be used in each direction without adding 
storage. 

No matter which way the spatial derivatives are computed, it is important 
to note that computing the Chebyshev derivative approximations requires only 
mesh point values. Derivatives at the end points require only points interior 
to the mesh so no extra procedure is required to compute derivatives at 
boundaries. 

Once the spatial derivatives are approximated, what results is a system of 
ordinary differential equations in time for the variation of the solution at 
each collocation point (Method of Lines). Because the differentiation matrix 
is full, explicit methods are typically used to integrate the semi-discrete 
equations. In this paper, all time integrations will be performed with a 
fourth-order Runge-Kutta method. 

The advantage of using this spectral method to solve (1) is that for 

GO 

solutions which are C (D), the accuracy is better than any polynomial order 
(Canuto and Quarteroni, [4]). This is usually called "spectral accuracy" and 
asymptotic behavior can be observed if there are enough grid points to 



-4- 


adequately resolve the solution. It is thus possible to compute to a given 
spatial accuracy with fewer grid points than required by typical low-order 
finite difference approximations. 

Balancing the high accuracy of the spectral method, however, are some 
major disadvantages of the typical Chebyshev collocation approach: 

(1) It may not be easy or even possible to map D D' globally. 

(2) The collocation point distribution is global and predetermined. Local 
refinement of the mesh is not possible. 

(3) The points are concentrated near the boundaries where they are 
typically not needed for hyperbolic problems. 

(4) If explicit time integration is used the time step restriction In one 
dimension is proportional to 1/N 2 . 

(5) For complete flexibility in the number of mesh points which can be 
used, the derivatives cost of 0(N 2 ) in each direction. 

These problems can be reduced significantly by breaking up the region D 
into several subdomains each of which has its own Chebyshev grid. With a 

stable and efficient method for computing the interfaces, the advantages of 
such an approach would be: 

(1) Complicated geometries can be covered. 

(2) Points can be distributed with some flexibility; local refinement is 
possible. 

one dimension, with N points and K subdomains, the time step 

2 

restriction increases to At « K/N . 

(4) Derivative evaluation work with matrix multiplication decreases to 
2 

K(N/K) or 1/K that of a single grid. 
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The idea of breaking up the computational domain into subdomains each with 
a different grid is not new. For finite difference methods this is a 
currently popular approach (e.g., (5]). For spectral methods, however, 
previous applications have been limited to elliptic and parabolic problems. 
Orszag [6] first applied such a technique to solve elliptic problems. He 
enforced continuity of the function and its first derivative as the interface 
condition. Metivet and Morchoisne [7] and later, .Morchoisne [8] computed 
multidomain solutions to the Navier-Stokes equations. Recently, Patera [9] 
and Korczak and Patera [10] have been using a spectral element method to solve 
the incompressible Navier-Stokes equations. Their method is very similar to 
the p finite-element methods developed by Babuska (see [10]) but uses 
Chebyshev interpolants. The treatment of the convective terms, however, does 
not lend itself to purely convective problems. For these problems, we 
describe the method below. 


2. MULTIDOMAIN APPROACH 

In this paper, we will break up the physical domain, D, into K 
subdomains which do not overlap except for the common boundary points. 

Figure 1 shows a rectangular two-dimensional example of the situation with 
four subdomains. Each of the are mapped onto a square [ — l,l]x[— 1,1] . 

Spatial approximations at interior points of each subdomain are computed in 
the usual way. Across an interface, however, there are two values of the 
normal derivative. For example, at the y coordinate line interface between 
Dj and D 2 in Figure 1, derivative approximations are available from the 


left and from the right. The problem Is to choose properly information from 



the right and the left to give a stable and consistent approximation to the 
differential equation at the interface. 

Before discussing a multidomain method for the boundary value problem (1), 
we will first examine the one-dimensional case. In one dimension, we seek 
interface algorithms of the semidiscrete form 




R 3Q R 
3x 


0 


(7) 


where Q denotes the value of Q at an interface and the derivatives 
superscripted with L and R denote the two spectral approximations computed 
in the left and right, respectively. For consistency, we require that 

A L + A R = A (8) 


and for efficiency we want A R and A R to be computed with little more work 
than is required for the computation of A itself. 

To generate the coefficient matrices, consider first the linear scalar 
hyperbolic equation 

U t + Xu x = 0 X > 0. (9) 


Because the equation is hyperbolic, it is clear that the common interface 
point should depend only on information propagated from the left. Thus, the 
approximation should be 

3u^ 3u^ 

3t A 3x 


0 . 


( 10 ) 
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This is, of course, just upwind differencing at the interface and is 
equivalent to the way Gottlieb and Orszag [11] handled a tau approximation to 
equation (9). To simplify the computational logic to include cases where the 
coefficient, X, is of either sign, the approximation (10) can be written as 

ff“ + V 2 (X + | X | ) |ij- + V 2 (X - | X | ) |^- - 0. (11) 

If we now consider that this equation is a single component of a 
diagonalized system, where the diagonal matrix 

<= P -1 AP, 

we can write the system as 

|^- + V 2 (A + |A| ) |5- + V 2 (A - |A| ) |2- =0 (12) 

where | A | = P | A | P Formally, this is nothing more than the method of 

characteristics in one dimension. 

We now propose to avoid the computation of the matrix absolute value by 
approximating it with a diagonal matrix 

] A j - PX* IP -1 = X* I (13) 

* 1 1 

where X is chosen to lie between the largest and smallest elements of |A|. 

The boundary scheme is now of the form of Eq. (7) with 
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A L - V 2 (A + X* I) A R - V 2 (A - X* I). (1A) 

This choice of coefficient matrices always has proper upwind dominance on 
characteristic variables, but includes some downwind influence. To 
see this, re-diagonalize the system (7) and use u as the n C ^ component of 
the diagonalized system. Then the approximation to the method of 
characteristics causes the characteristic variables at the interface to be 
approximated by 


3u 

3t 






0. 


(15) 


In fact, this can be viewed as the purely upwind scheme with an error term: 

For the X >0 case. 

n * 


3u^~ 3u^ 

3 1 n 3x 




3x 


(16) 


Thus, we have the spectrally accurate upwind approximation with an error 
term proportional to the difference of the right and left spectral 
derivatives. If the solution has the necessary smoothness, this difference 
should also decay spectrally and spectral accuracy of the approximation should 
be retained. 

We will study the stability of the multidomain method with the interface 
approximation (1A) numerically. An analytic study of stability is not 
possible at this time. Stability theory for Chebyshev approximations to 
hyperbolic initial-boundary value problems is not advanced enough to analyze 
an approximation which introduces some downwind influence at the interface. 



-9- 


We consider the two-domain approximation of the scalar equation (9) with 
the interface approximation (12) with X ■ 1. The line segment [-2,2] is 
divided equally into two domains of [-2,0] on the left and [0,2] on the right. 
The semidiscrete approximation can be written as a system of ordinary 
equations with the two-domain coefficient matrix 


D 



0 



(17) 


L R 

where D and D are the single domain differentiation matrices for the 
left and the right, modified to include the interface approximation. For this 
system to be time stable, that is, the solution does not grow unboundedly as 
t the eigenvalues of the coefficient matrix must have negative real 

parts . 

* 

Figure 2 shows how the eigenvalues change as X varies when 6 points are 
* 

used. The case of X = 0 corresponds to simple averaging and is clearly not 

* 

time stable. Choosing X >0 large enough moves the eigenvalues into the 
left half of the complex plane and the resulting approximation is time 

rk 

stable. The case of X = 1 is the purely upwind case and the eigenvalues 

* 

decouple into two single-domain patterns. If X is chosen equal to, or 

larger than, the wave speed, X n , the approximation has the effect of adding a 

purely dissipative term to the equation and two purely real eigenvalues are 
* 

created. If X is very much larger than X n , however, the eigenvalues 

migrate to the right of the imaginary axis. The range of X's for which the 
approximation is stable decreases as the disparity in the number of points 
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becomes larger; for very stiff systems, it may be necessary to use |a| 

* 

instead of X at the interface. 

It is interesting to note that the reverse situation, where there is more 

resolution on the upstream side of the interface, does not show this behavior 

* * 

and is stable for all X _> 0. For systems, this means that X should be 

chosen to be only slightly larger than the smallest eigenvalue representing a 

characteristic moving from the coarse to the fine grid. For systems, this 
* 

means that X should be chosen to be only slightly larger than the smallest 
eigenvalue representing a characteristic moving from the coarse to the fine 
£rid. We note, however, that the examples on which the scheme has been tested 
show that the approximation is robust over a wide range of choices of X*. 

In two dimensions, the upwind weighted approximation is used in the 
direction perpendicular to the interface. Returning to Figure 1, along x 
coordinate lines, the y derivatives are continuous across the interfaces 
except at corners. At points not on the corners, then, we propose using 


3Q 

3t 


I 


+ 



+ A 


R 3Q R 
3x 


+ 



(18) 


L# R 

where A and A are defined as above. Along x coordinate interfaces, 


|2i + A >a! + b t |q!i + b b |q! 

3t 3x 3y 3y 


0 


(19) 


where B T = l / 2 (B + V* I) and B B = V 2 ( B - y* I) and y* is an approxima- 
tion to the eigenvalues of B. At corners ,- the weighted approximations are 
used in both directions. 
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3. NUMERICAL EXAMPLES 

Numerical experiments on four model problems in one and two dimensions 
will be presented. The models include the scalar one-dimensional hyperbolic 
initial boundary value problem for a travelling Gaussian pulse, a linear 
system in one dimension, quasi-one-dimensional flow in a converging-diverging 
nozzle, and the transonic Ringleb problem. The Ringleb flow models the smooth 
nonlinear transonic flow in a curved duct and has an exact solution to which 
to compare. 

A. Solution of a Linear Scalar Problem 

The solution to the linear scalar problem 

+ 2 — = 0 xe [-2,2], t > 0 (20) 

u(x,0) = exp(-(x - Xq) 2 /0.3) x € [-2,2]- 

u(-2,t) - exp(-(x - t - Xq)- 2 /0.3) t > 0 

* 

can be used to examine the effects of varying X in the spatial 

approximation described in Eq. (15). The time integration for this and all 

following examples was a fourth-order Runge-Kutta technique. For this and the 

next model problem the time step was chosen so that the temporal errors were 

on the order of 10 - ^. The main questions to be answered here are the effect 
★ 

of the A * 2 on the accuracy of the solution and if reflections are a 

problem at - the interface. Figure 3 shows the computed (circles) and exact 

(line) solutions for the pulse after it has propagated through the interface 

* 

at x-= 0 for two distributions of the mesh points and X 


- 6 . 
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The interface approximation Eq. (15) degrades the accuracy of the solution 

when compared to the purely characteristic interface, X •* 2, if equal 

resolution is not provided in each subdomain. In no case, however, is the 

global 1.2 error larger than the global error for the characteristic inter— 

* 

face. Furthermore, if X remains fixed and the total number of points is 

increased, the error decay remains spectral. Figure 4 shows the pointwise 

errors of the solution to Eq. (20) for the situations represented in Figure 3 
* 

as X is increased beyond the characteristic value of 2. The situation is 
worse when more resolution is used upstream of the interface because the 
approximation includes more and more downwind influence as X* is 
increased. In a practical computation, the effect of the boundary 
approximation would not be important if the solution were equally resolved in 
all subdomains. 

Reflections at the interface are not visible in Figure 3 even though there 
is a factor of two difference in the number of collocation points. Gottlieb 
and Orszag [11] also noticed this for a tau approximation to the scalar wave 
equation. This is typical for the spectral approximations; examples with up 
to a factor of three and four in the ratio of the number of mesh points have 
not shown spurious reflections off of the interface. 

B. A Linear System Example 

The accuracy of the interface approximation will now be demonstrated with 
the 2x2 linear system 

x £ [-2,2], t > 0. (21) 

x 
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The coefficent matrix has eigenvalues +3 and -1 so the system has 

information which propagates in both directions and with different speeds 

across the interface at x «= 0. The initial and boundary conditions were 

chosen so that the characteristic variables were the Gaussian pulses used in 

* 

the scalar problem, Eq. (20). The coefficient X for this case was chosen 
to be the maximum eigenvalue, X = +3. Figure 5 shows the results for the two 
components of this system at a time when the characteristic pulses have 
crossed the interface. In Figure 5a there are twice as many points to the 
left of the interface as to the right and this is reversed for Figure 5b. The 
symbols represent the computed solutions and the solid lines represent the 
exact solutions. 

A study of discrete L£ errors for the system computations is shown in 
Tables I through III. Clearly, the error is spectral for all three 
situations. In fact, for an equal number of mesh points on either side of the 
interface, the error decay is exponential. For the problem of propagating 

pulses, where the features needing higher resolution are continually moving, 
it is not surprising that the best errors are obtained when there are an equal 
number of mesh points on both sides of the interface. 


C. Quasi-One-dimensional Nozzle Flow 

One potential point of concern in using the interface approximation given 
by Eq. (14) regards the stability of cases where one of the eigenvalues of the 
coefficient matrix is much larger than any other. Such a situation occurs at 
sonic points in an ideal gas flow where one of the characteristic speeds. 


actually vanishes 
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To test this situation the nonlinear problem of steady gas flow in a 
quasi-one-dimensional converging-diverging nozzle was solved with the 
multidomain method where an interface was placed at the sonic point. The 
quasilinear form of the Euler gas dynamics equations for time-dependent flow 
in a quasi— one— dimensional nozzle without shocks can be written as 


■p" 


u y‘ 


'P' 


yuA (x)/A(x) 


+ 




s 

X 

u 

t 

2, 

a /y u 


u 


0 

- 


. 


. J 

X 



where P is the logarithm of the pressure, u is the gas velocity, y is the 
ratio of specific heats, and a is the sound speed. The coefficient matrix 
has eigenvalues of u + a and u — a so that one of them is zero at a sonic 
point. The steady flow is found as the large time limit of the unsteady flow 
described by (22). 

The nozzle area is given by A(x) = x/2 + 1/x so the throat occurs at 
x = ^2. For the cases run, a subsonic inflow boundary was placed at x = 0.2 
and characteristic boundary conditions were used. After the gas accelerates 
through the sonic value at the throat, it leaves the nozzle supersonically so 
no boundary conditions are applied at the outflow. 

For the gas dynamics calculations in one dimension, X = ^2(l u+a l + |u-a|) 
was chosen since this corresponds to the diagonal elements of the absolute 
value of the coefficient matrix. Although the problem was solved for domain 
interfaces in both the subsonic and supersonic portions of the nozzle, only 
results for a single interface at the sonic point will be shown here. (The 
two-dimensional example below will include a variety of interface placements.) 
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Figure 6 shows the steady pressure in the nozzle computed with two domains 
and twice as many mesh points on the right as on the left. Our tests on a 
variety of grids have not shown any stability difficulties in computing steady 
flows when placing the interface at a sonic point. 


D. Two-Dimensional Transonic Flow 

A more complicated problem is the two-dimensional transonic Ringleb 
flow. This problem allows us to study the computational efficiency of the 
multidomain solution algorithm as outlined in the Introduction. Kopriva, et 
al., [12] used this problem for a comparison of the performance of the 
spectral method with a second-order finite-difference method. In this section 
we will compare the multidomain spectral method with the single domain 
spectral method. 

The Ringleb flow is a simple example of a two-dimensional transonic flow 
for which there is an exact solution. (See, for example, Courant and 
Friedrichs [13].) The streamlines of the physical space solution appear at 
large distances as parabolas which are determined from a special hodograph 
solution of the potential equation for steady irrotational isentropic flow. 
By choosing two streamlines to represent solid walls, this problem models a 
steady transonic flow in a duct. Figure 7 shows the Mach contours of one such 
duct flow. 

Again we will look for the large time solution of the unsteady gas 
dynamics equations, this time in two dimensions. The problem in the curved 
duct shown in Figure 7 is mapped onto a rectangle in the stream function- 
potential coordinate system derived from the exact solution. In this 
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coordinate system, the unsteady equations can be written as 


Q t - -R (23) 


where R is the steady state residual 


R - A % + B V (24) 


Since the solution is irrotational, the solution vector is chosen to be 


Q - [P u v] T (25) 


and the coefficient matrices are 


u 

♦ 

X 

<p 

y 

0 


V 

i|» 

X 

♦ 

y 

o’" 

a 2 t x /y 

U 

0 

0 

B = 

a 2 $ x /y 

V 

0 

0 

a 2 4> y /Y 

0 

u 

0 


2 

a i|< /Y 

0 

V 

0 

0 

0 

0 

u 


0 

0 

0 

V 

— 





— 





As before, P represents the logarithm of the pressure and (u,v) represent 

the velocity components in the Cartesian x and y directions, 

respectively. The matrix coefficients are computed from the mapping derived 

from the exact solution and the contravariant velocity components are 

U = u<J) + vd> and V = uib + vib . 

Y x T y T x y 
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The physical boundary conditions for this problem represent subsonic 
inflow at the entrance of the duct (at the lower left of Figure 7), supersonic 
outflow at the exit, and the sides are treated as impermeable boundaries 
(walls). So that the initial boundary value problem is well-posed the 
boundary conditions must be chosen carefully. See Kopriva, et al., [12] for 
details of the procedure which follows. For the subsonic inflow, we can 

specify only two quantities and have chosen the total enthalpy and the angle 
of the flow (so V - 0). The quantities P and U are computed from two 
conditions: The first is a compatibility equation derived from the pressure 

equation and the normal momentum equation. The second comes from 
differentiating the enthalpy equation in time. From U and the condition 
V = 0, the Cartesian velocities u, v can be computed. At the outflow, no 
boundary conditions are needed. Finally, at the walls the normal velocity, U, 
must vanish. The vector Q is computed by solving the tangential momentum 
equation for V and a compatibility equation which combines the normal 

momentum and pressure equations for P. 

The system of equations (22) were discretized ' as described above, and 
fourth-order Runge-Kutta was used for the time integration. For a single 
domain, the Chebyshev spectral grid for the Ringleb problem with 16 streamwise 
and 8 normal mesh intervals is shown in Figure 8. It is clear that the 
spectral method strongly concentrates the grid points near the walls. The 
largest gradients, however, occur in the streamwise direction near the sonic 
line (as can be seen in Figs. 7 and 9) where the streamwise mesh distribution 
is coarsest. These two factors contribute to the fact that the time 
integration step is very small and that accuracy is degraded by the lack of 
resolution where it is needed. 



-18- 


A multidomain grid distribution for which performance will be compared to 
the single domain method is shown in Figure 10. Six domains now cover the 
duct and the same number of mesh intervals as for the single domain case are 
used. The divisions were chosen to demonstrate the kinds of situations which 
the multidomain method should be able to handle. Three divisions with 
6 + 5+5 mesh intervals are in the streamwise direction and two are in the 
normal direction. With this choice, two points occur where the corners of 
four domains come together. The first domain boundary in the streamwise 
direction was chosen to appear in a subsonic region of the duct. The second 
domain boundary in the streamwise direction was chosen to intersect the sonic 
line. By dividing the normal direction into two domains, the effective mesh 
spacing near the walls is doubled. Finally, note that by comparing Figure 10 
to Figure 7 the sonic line also intersects the domain interface in the normal 
direction. 

To allow comparison. Figure 11 shows the Mach number contours for both the 
single domain and the multidomain solutions. Note particularly that the sonic 
line remains smooth through the domain interfaces. Table IV summarizes the 
performance of the single domain spectral method compared with this particular 
choice of grid. First, note that even with this distribution of domains, the 
maximum error in the pressure for the multidomain computation has not been 
degraded from the single grid one. In fact, the error is five percent better. 

The real advantage that the splitting has had for this case, however, is 
that the multidomain solution relaxes more quickly to steady state for a given 
number of intervals and accuracy. Figure 12 compares the rate at which the 
discrete 1^ norm of the residual of the pressure decays for the single and 
multidomain cases. The results are also summarized in Table IV. From the 
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trend of the graph, it should take over 21/2 times as many iterations for the 

single grid residual to decay to that of the single grid residual. This is a 

direct result of the fact that larger time steps can be used for the multi- 

* 

domain case. The choice of X also affects the convergence rate: larger 

values up to the stability limit give faster convergence to steady state. 

The advantage of a k-domain derivative computation requiring 1/k the 
amount of work as a single domain computation does not show up in this 
example. In fact, as Table IV indicates, the average time per iteration (time 
step) requires the same amount of time at 0.5 sec. on the Langley Cyber 855. 
This is due to the fact that there is overhead in computing the interface 
approximation. Doubling the number of points in each direction with the same 
domain distribution decreases the time per iteration for the multidomain 
computation to 70% of the single domain cost. Though no attempt was made to 
compute the interface conditions efficiently, the number of points inside each 
domain will have to be large compared to the number of domains for the 
efficiency gained by being able to use fewer points in computing derivatives 
to become important. 

The final advantage of a multidomain method which was listed in the 
Introduction is that flexibility in the choice of grid point distribution is 
now possible. A series of calculations were made with the duct being divided 
into two domain intervals in each direction. As with Figure 10, the direction 
across the duct was divided in half and the same number of mesh points was 
used. In the streamwise direction, however, only one domain boundary was 
inserted. This boundary was inserted in several places along the duct with 
different numbers of points on either side. 
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Results of some of the computations are summarized in Table V. The 
division is reported in terms of the fraction of the total variation of the 
velocity potential along the length of the duct. The first entry in the list 
places the division approximately near the bend of the duct where the 
gradients of the solution are the highest. It is clear that with a proper 
choice of grid it is possible to obtain better accuracy with the multidomain 
distribution of a given number of grid points than with a single grid. For 
the best case computed here, the error is about 2 1/2 times better for the 
multidomain calculation. 

The problem of how to properly distribute points and subdomains in general 
is a major one and is beyond the scope of this paper. If they are poorly 
placed the error can be worse than the single domain error (see Table V). For 
now, it is not known how to obtain the optimal point and subdomain 
distribution. Rather, some knowledge of the behavior of the solution must be 
used as a guide. 


CONCLUSIONS 

We have described a simple approximation which allows a multidomain 
spectral solution of quasilinear hyperbolic equations. Numerical examples of 
linear equation models and ideal gas flow show that the method gives 
advantages in both accuracy and efficiency over using a single domain. 
Dividing up a computational domain into several subdomains gives the 
possibility of local refinement and allows some flexibility in the 
distribution of mesh points. It is possible to obtain better accuracy by 
doing so. Also, with multiple domains it is possible to take larger time 
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steps than with a single domain. This increases the efficiency for using time 

relaxation to acheive steady state solutions. 

The use of a multidomain technique is also appropriate if discontinuities 
are fitted as boundaries. When shocks occur within a flow, subdomains would 
be arranged so that each shock lies on a subdomain boundary. In smooth parts 
of the solution, the technique described here would be used. Along shock 
interfaces, a shock fitting algorithm like that described in reference [1] can 
be used (Kopriva and Hussaini, to be published). 

The theoretical issues which remain are many. Some theory for the range 
of values which X can take for the method to be stable must be found. 
However, choosing X to be the average of the largest and smallest 
eigenvalues of the coefficient matrix has always worked. Finally, like the 
problems associated with the p — version of the finite-element method, the 
choice of domain and point distribution for a given number of points is an 
open issue. 
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TABLE I. l >2 errors for the solutions to Eq. (20) with equal 
number of points on each side of the interface. 


N 

Error in u 

Error in v 

8 

1.57 x 10" 2 

1.49 x 10 -2 

16 

4.15 x 10" 6 

4.86 x 10 -6 

32 

1.91 x 10~ 9 

1.91 x 10" 9 


TABLE II. L2 errors for the solutions to Eq. (20) with more 
points to the right of the interface. 


N L’ N R 

Error in u 

Error in v 


-2 

-2 

8, 16 

1.22 x 10 

1.05 x 10 


-4 

-4 

12, 24 

2.45 x 10 

2.33 x 10 


-6 

-6 

16, 32 

3.93 x 10 

3.93 x 10 


TABLE III. L2 errors for the solutions to Eq. (20) with 
more points to the left of the interface. 


N L’ N R 

Error in u 

Error in v 


-3 

-2 

16, 8 

9.80 x 10 

1.04 x 10 


-4 

-4 

24, 12 

3.48 x 10 

2.88 x 10 


-f, 

-6 

32, 16 

> 

O 

X 

• 

fH 

2.30 x 10 
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TABLE IV. Performance comparison for single and multidomain spectral 
computations. 

Grids: Single Domain (SD) 17 x 9 points 

Multidomain (MD) (7 + 6 + 6) x (5+5) points 
(separated by domain) 

Maximum Error 
SD 
MD 

Number of Steps to Reduce Residual Three Orders of Magnitude 
SD > 1500 

MD 780 

Average Spectral Radius 

SD 0.9964 

MD 0.9942 

Average Time per Iteration 

SD 0.50 sec. 

MD 0.50 sec. 


1.85 x 10 
1.74 x 10' 


TABLE V. Effect of streamwise mesh distribution 
on Ringleb calculation. 


Grid 

Division 

Maximum Error 

8 + 8 

0.45 + 0.55 

7.8 x 10" 4 

8 + 8 

0.50 + 0.50 

9.3 x 10" 4 

16(SD) 

— 

1.9 x 10~ 3 

10 + 6 

0.34 + 0.66 

1.2 x 10“ 2 
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x — > 


FIG. I. Diagram of the two-dimensional subdomain structure used to divide a 
computational domain. 













[012 

X 

3a. Solution of the scalar pulse problem Eq. (19) computed on two 
domains shown after the pulse has travelled from the left through 
the Interface at x ° 0. Computations are for 22 points left and 11 
points right of the interface. The exact solution is the solid 


line; computed solutions are the circles. 


0 1 2 

X 

Solution of the scalar pulse problem Eq. (19) computed, on two 
domains shown after the pulse has travelled from the left through 
the interface at. . x = 0. Computations are for 11 points left and 22 
points right. The exact solution is the solid line; computed 
solutions are the circles. 
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FIG. 4a. Pointwise errors as X varies for the situation in Figure 3a. 


- 2.0 - 1.5 - 1.0 --5 0 .5 1.0 1.5 2.0 

X 

* 

FIG. 4b. Pointwise errors as X varies for the situation in Figure 3b. 



-2 -1 01 2 

X 

FIG. 5a. Graph of the two solutions u (circles) and v (squares) of the 
linear system Eq. (20) with 22 points on the left and 11 points on 
the right. The exact solutions are represented by the solid line. 
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-2 -1 0 1 2 


X 

FIG. 5b. Graph of the two solutions u (circles) and v (squares) of the 


linear system Eq. (20) with 11 and 22 points on the left and the 
right. The exact solutions are represented by the solid line. 
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FIG. 8. Single domain Chebyshev grid for the Ringleb problem. 
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FIG. 10. Multidomain grid with six subdomains for the Ringleb problem. 
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FIG. 12. Comparison of residual decay for single domain and multidomain 
solutions to the Ringleb problem. 
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